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A  simulation  algorithm  is  presented  for  multi-agent  hybrid  systems — systems  consisting  of  many 
sets  of  nonsmooth  differential  equations — such  as  systems  involving  multiple  rigid  bodies,  vehicles, 
or  airplanes.  The  differential  equations  are  partitioned  into  coupled  subsystems,  called  “agents”; 
and  the  conditions  which  trigger  the  discontinuities  in  the  derivatives,  called  “events”,  may  depend 
on  the  global  state  vector.  Such  systems  normally  require  significant  computational  resources  to 
simulate  because  a  global  time  step  is  used  to  ensure  the  discontinuity  is  properly  handled.  When 
the  number  of  systems  is  large,  forcing  all  system  to  be  simulated  at  the  same  rate  creates  a 
computational  bottleneck,  dramatically  decreasing  efficiency.  By  using  a  control  systems  approach 
for  selecting  integration  step  sizes,  we  avoid  using  a  global  time  step.  Each  subsystem  can  be 
simulated  asynchronously  when  the  state  is  away  from  the  event.  As  the  state  approaches  the 
event,  the  simulation  is  able  to  synchronize  each  of  the  local  time  clocks  in  such  a  way  that  the 
discontinuities  are  properly  handled  without  the  need  for  “roll  back”.  The  algorithm’s  operation  and 
utility  is  demonstrated  on  an  example  problem  inspired  by  autonomous  highway  vehicles.  Using 
a  combination  of  stochastic  modelling  and  numerical  experiments  we  show  that  the  algorithm 
requires  significantly  less  computation  time  when  compared  with  traditional  simulation  techniques 
for  such  problems,  and  scales  more  favorably  with  problem  size. 
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1.  INTRODUCTION 

Many  engineering  systems  are  best  described  by  sets  of  ordinary  differential 
equations  (ODEs)  with  discontinuous  right-hand  sides.  Examples  include  phase 
transitions,  contact  mechanics,  or  the  dynamics  of  physical  systems  controlled 
by  digital  computers.  Such  systems  arise  in  many  contexts  and  are  often  re¬ 
ferred  to  as  hybrid  systems,  switched  systems,  or  nonsmooth  systems.  They  can 
also  be  though  of  as  discrete  event  systems  augmented  with  differential  equa¬ 
tions.  See  Branicky  [1995]  for  a  discussion  of  various  modeling  paradigms.  The 
most  basic  form  of  such  a  systems  is 

[  fa(x,t),  g(x)  <  0 

X  =  {  (1) 

I  fb(x,  t),  g{x )  >  0. 


Each  of  the  functions  fa  and  ft,  describing  the  derivative  represent  distinct 
modes  of  operation  of  the  hybrid  system.  We  will  refer  to  each  function  simply 
as  a  mode. 

Changes  from  one  mode  of  operation  to  another,  called  mode  switches,  are 
caused  by  state  events  or  simply  events.  The  occurrence  of  these  events  is  trig¬ 
gered  by  the  zeros  of  an  event  function  g(x(t))  (also  called  a  guard  or  discon¬ 
tinuity  function  in  the  literature).  When  such  an  event  occurs  the  first  deriva¬ 
tive  of  the  solution  trajectory  becomes  discontinuous.  Of  particular  interest  is 
the  fact  that  embedded  or  software  based  control  systems  can  be  modeled  in 
such  a  fashion  (see,  e.g.,  Maler  and  Pnueli  [2003]  and  other  in  that  series). 
In  such  a  case,  there  can  be  many  modes  and  their  connectivity  may  be  quite 
complex. 

Multi-agent  hybrid  systems  are  groups  of  individual  interacting  hybrid  sys¬ 
tems.  More  precisely,  we  define  multi-agent  hybrid  systems  as  collections  of 
individual  hybrid  systems,  each  called  agents,  in  which  the  there  is  no  agent- 
to-agent  coupling  in  the  right  hand  sides  of  the  differential  equations.  However, 
the  event  functions  which  trigger  mode  transitions  may  depend  on  the  global 
state  vector.  This  concept  is  best  illustrated  through  an  example.  Consider  the 
simplified  kinematics  of  two  automated  highway  vehicles: 
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where  the  superscripts  1  and  2  are  used  to  distinguish  between  the  state  (posi¬ 
tion  and  orientation  in  the  plane)  of  the  first  and  second  cars  (i.e.,  agents  1  and 
2),  x1  =  [ z 1  y1  d1]  andx2  =  [z2  y2  02].  The  functions  ux(t),  v2(t),  co\t)  and  w2(t) 
are  the  forward  and  turning  velocities  of  the  two  cars  respectively.  These  time- 
dependent  (perhaps  implicitly)  functions  are  considered  inputs  to  the  system 
and  are  supplied  by  the  automatic  control  system.  The  control  system  may  have 
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a  library  of  possible  functions  for  the  velocity  and  turning  rate,  which  it  selects 
among  based  on  current  operating  conditions.  For  example,  under  nominal  op¬ 
eration  the  vehicles  are  to  drive  in  a  straight  line;  however,  if  the  two  vehicles 
come  within  R  feet  of  each  other,  an  emergency  collision  avoidance  controller 
would  be  activated  and  both  the  cars  would  steer  away  from  one  another.  In 
such  a  situation  perhaps: 


!  _  Jo,  gix1,*2)  >  0 

I  1,  gix1,  x2)  <  0 

gix1,  x2)  =  (x1  -  x2)2  +  (y 1  -  y2)2  -  R2  <  0. 


(4) 

(5) 


Obviously  much  more  sophisticated  examples  and  control  laws  could  be  used, 
however,  the  representative  features  of  this  example  are  that  the  right-hand 
sides  of  the  differential  equations  for  an  individual  car’s  dynamics  do  not  depend 
on  the  state  of  the  other  car;  yet  the  condition  for  switching  control  laws  depends 
on  the  states  of  several  agents. 

Examples  of  such  systems  abound.  Automated  highway  systems  [Deshpande 
and  Semenzato  1995],  free-flight  air  traffic  control  [Tomlin  et  al.  1998],  coop¬ 
erative  mobile  robotic  systems  [Fierro  et  al.  2002],  cellular  biology  [Alur  et  al. 
2001],  and  multibody  systems  simulated  for  graphics  applications  all  can  be 
modeled  this  way — the  dynamics  of  each  vehicle,  plane,  or  robot  are  decoupled; 
however,  certain  critical  events  which  are  relevant  to  the  simulation  (e.g.,  col¬ 
lisions)  depend  on  the  states  of  pairs  of  agents.  Also,  note  that  the  partitioning 
of  the  agents  can  change  throughout  the  course  of  the  simulation.  Two  bodies 
can  collide  and  stick  together  for  a  brief  duration,  in  which  case  their  ODE’s 
are  coupled  and  they  are  a  single  agent.  Likewise,  the  collision  avoidance  mode 
for  the  automated  vehicle  example  may  make  use  of  a  closed  loop  feedback  law, 
in  which  case  the  velocity  for  car  1  is  determined  based  on  the  position  of  car  2. 
Coupling  of  the  differential  equations,  either  through  software  or  physical  in¬ 
teraction,  makes  it  necessary  to  consider  the  two  systems  as  subsystems  of  a 
single  agent  for  the  purposes  of  simulation. 

Traditionally,  the  need  to  properly  detect  and  manage  such  state  events  has 
suggested  the  use  of  a  single  global  integration  step  size,  which  is  necessarily 
the  minimum  acceptable  step  size  across  all  agents.  Because  the  differential 
equations  of  each  agent  are  decoupled,  this  is  not  required  and  often  leads 
to  severe  loss  of  simulation  efficiency,  especially  in  applications  with  a  large 
number  of  agents.  This  phenomena  is  explained  in  Section  2. 

This  article  introduces  an  alternate  simulation  algorithm  that  allows  each 
agent  to  use  a  different  step  size,  resulting  in  an  asynchronous  simulation. 
The  effect  is  that  the  average  step  size  is  much  larger  resulting  in  significant 
decreases  in  computation  time.  In  Section  3,  we  review  work  related  to  asyn¬ 
chronous  simulation.  In  Section  4,  our  methodology  for  selecting  integration 
step  sizes,  while  properly  handling  state  events,  is  presented  along  with  re¬ 
sults  on  the  need  for  “roll  back”.  The  overall  algorithm  is  introduced  in  Section  5 
and  illustrated  via  example  in  Section  6.  Section  7  quantifies  the  algorithms 
advantages  over  the  traditional  approach  in  terms  of  computational  cost.  The 
role  of  the  step  size  selection  scheme  is  modelled  statistically  and  the  effect  of 
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the  additional  overhead  of  the  new  algorithm  is  examined.  In  Section  8,  the 
performance  improvement  is  determined  experimentally. 

2.  BACKGROUND  AND  CHALLENGES 

In  Cellier  [1979],  it  was  shown  that  the  proper  way  to  simulate  any  hybrid 
system  is  to  numerically  integrate  all  of  the  differential  equations  until  t*  — 
the  first  time  at  which  g(x(t))  =  0.  At  this  point,  the  numerical  integration  is 
stopped,  any  applicable  mode  switches  are  activated  and  the  integration  may 
be  restarted,  using  x(t*)  as  the  new  initial  condition.  This  technique  is  referred 
to  as  discontinuity  locking  and  is  widely  accepted  as  the  standard  hybrid 
system  simulation  methodology.  This  requirement  of  stopping  the  integration 
precisely  when  g(x(t*))  =  0  gives  rise  to  the  event  detection  problem.  Since 
numerical  integrations  are  performed  in  discrete  time,  T  e  {t\,  4,  ■  ■  . ,  4,  •  ■  ■} 
with  step  sizes  hi,  =  4  —  4- l,  it  is  difficult  to  find  t*  exactly.  Much  work  has 
been  done  on  the  problem  and  most  reliable  approaches  use  interpolants  to 
approximate  the  state  between  steps  and  then  check  these  interpolants  to  find 
the  time  of  zero  crossings  of  g(x )  (e.g.,  Shampine  et  al.  [1991],  Park  and  Barton 
[1996],  Esposito  et  al.  [2001b],  and  Bahl  and  Linninger  [2001]). 

It  is  well  known  that,  when  simulating  hybrid  systems,  a  failure  to  detect  an 
event  can  have  disastrous  results  on  the  global  solution  due  to  the  nonsmooth 
nature  of  the  problem  [Branicky  1995].  Works  detailing  requirements  for  hy¬ 
brid  simulators  list  accurate  event  detection  as  a  primary  concern  [Mosterman 
1999].  Event  detection  in  hybrid  systems  is,  in  itself,  a  very  difficult  prob¬ 
lem  [Shampine  et  al.  1991]. 

A  second  important  issue  in  numerical  integration  is  managing  the  tradeoff 
between  efficiency  and  the  desired  accuracy,  on-the-fly,  by  varying  the  step  size. 
Large  step  sizes  result  in  more  efficient  simulations  but  decrease  the  accuracy 
and  stability  of  the  results.  Small  step  sizes  are  required  to  maintain  strin¬ 
gent  accuracy  requirements  when  the  solution  to  the  differential  equation  is 
ill-behaved,  however  the  computations  become  extremely  expensive  due  to  the 
large  number  of  steps  needed  to  complete  the  simulation.  Automatic  step  size 
selection  schemes  attempt  to  optimize  this  tradeoff  by  selecting  the  largest 
possible  step  size  consistent  three  criteria:  (1)  maintaining  the  estimated  trun¬ 
cation  error  within  a  desired  range;  (2)  maintaining  numerical  stability;  and 
(3)  accurate  event  detection. 

Traditionally,  all  hybrid  system  simulators  use  a  single  global  notion  of  time, 
meaning  the  integration  of  all  agents  is  synchronized  using  the  same  time 
discretization.  In  such  a  scheme,  at  each  iteration,  a  candidate  integration  step 
size,  h‘  for  agent  i ,  is  computed  for  each  agent.  This  candidate  represents  the 
maximum  step  which  will  result  in  an  acceptable  truncation  error.  The  only 
acceptable  choice  of  a  global  step  size-if  one  wishes  to  maintain  each  agent 
with  the  acceptable  range  of  truncation  error-is  the  minimum  among  all  the 
agent’s  candidate  step  sizes,  h  =  minf/i1, . . . ,  hN],  This  selection  means  some 
agents  will  be  simulated  with  an  unnaturally  small  step  size — resulting  in  a 
generally  inefficient  integration  scheme.  In  particular  as  the  number  of  agents 
increases,  h  tends  to  decrease  on  average. 
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Note  that  the  decision  to  use  a  global  step  size  is  not  motivated  by  the  differ¬ 
ential  equations.  Had  this  been  a  purely  smooth  system  of  decoupled  differential 
equations  (no  events),  the  optimal  efficiency  would  be  obtained  by  permitting 
each  agent  to  use  the  largest  acceptable  step  size.  This  would  result  in  an 
asynchronous  simulation  in  which  each  agent  proceeded  at  its  own  “fastest” 
acceptable  simulation  rate. 

The  motivation  for  using  a  traditional  global  choice  of  step  sizes  is  entirely 
attributed  to  the  coupled  nature  of  the  event  detection  problem.  If  agents  1 
and  2  are  simulated  with  different  step  sizes,  and  hence  different  time  meshes 
T 1  =  R* ,  f g  >  •  •  •  >  tl)  and  T2  =  {t2,  t2,  ■  ■  ■ ,  t2},  state  estimates  for  the  two  agents 
are  produced  at  different  points  in  time.  Deciding  if  an  event  has  occurred  now 
becomes  nontrivial  because  merely  evaluating  g(x 1,x2)  <  0  requires  x1  and 
x2  be  evaluated  at  the  same  time  instant — that  is,  gCxHfa), x2{tk))  is  mean¬ 
ingful  while  g(x1(tj),  x2(tk))  with  tj  ^  4  is  not.  Rather  than  address  this 
issue  traditional  simulators  simply  use  a  global  notion  of  time  to  simplify 
bookkeeping,  despite  the  decrease  in  performance.  Satisfying  these  two  com¬ 
peting  objectives,  increasing  efficiency  through  the  largest  possible  step  sizes 
while  ensuring  events  are  properly  detected,  is  the  central  challenge  addressed 
here. 


3.  RELATED  WORK 

Excellent  works  on  the  event  detection  problem  in  general  are  cited  above;  how¬ 
ever,  it  is  only  recently  that  specialized  modelling  languages  [Alur  et  al.  2000] 
and  simulation  environments  [Deshpande  and  Semenzato  1995]  have  become 
available  for  large  scale  multi-agent  hybrid  systems.  Relevant  work  includes 
two  tangentially  related  fields:  Multirate  numerical  integration  methods  and 
Distributed  Discrete  Event  System  Simulation  (DDESS). 

Given  a  group  of  coupled,  purely  smooth  differential  equations  (i.e.,  no 
events)  which  inherently  evolve  at  different  time  scales,  Multirate  integration 
techniques  attempt  to  efficiently  integrate  the  system  of  ODE’s  by  using  dif¬ 
ferent  integration  rates,  or  step  sizes  for  each  equation  (see  Gear  and  Wells 
[1984],  Engstler  and  Lubich  [1996],  or  Esposito  and  Kumar  [2001],  for  exam¬ 
ple).  By  using  the  appropriate  size  time  step  for  each  ODE  it  is  hoped  that  the 
overall  computational  effort  for  the  system  will  be  reduced.  Significant  chal¬ 
lenge  lies  in  attempting  to  simulate  two  coupled  differential  equations  using 
different  time  meshes.  Often  one  must  evaluate  the  derivative  of  the  so-called 
fast  variables  at  a  point  in  time  at  which  values  of  the  slow  variables  have  not 
been  computed.  This  is  typically  accomplished  by  constructing  some  appropri¬ 
ate  interpolant  for  the  slow  variables  so  that  they  may  be  approximated  at  off 
mesh  points,  when  needed. 

In  spirit,  the  goal  of  multirate  simulation  is  the  same  as  the  goal  of  this 
work:  reduced  computation  time  through  asynchronous  integration.  However, 
in  practice  the  challenges  are  quite  different.  The  primary  concern  of  multirate 
methods  is  accommodating  coupling  between  the  right-hand  sides  of  two  ODEs 
being  integrated  at  different  rates-a  nonexistent  issue  in  the  multiagent  prob¬ 
lem  because  all  coupled  ODE’s  in  our  scheme  are  simulated  at  the  same  rate. 
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Likewise  multirate  methods  do  not  address  the  problem  of  locating  discrete 
events — the  central  concern  in  this  article. 

A  second  more  closely  related  area  is  Distributed  Discrete  Event  System  Sim¬ 
ulation  (DDESS).  Given  a  set  of  loosely  decoupled  but  interacting  Discrete 
Event  Systems  (DES),  the  idea  is  to  make  the  best  use  of  any  available  compu¬ 
tational  resources  by  simulating  each  DES  on  a  different  computer  on  a  local 
network.  Each  DES  is  referred  to  as  a  “process”  and  is  simulated  in  isolation. 
When  a  DES  generates  a  message  which  is  relevant  to  another  DES,  a  message 
is  sent  over  the  network  to  the  computer  that  is  simulating  that  process.  Com¬ 
plications  arise  because  at  any  given  instant,  some  of  the  individual  process 
simulations  may  have  progressed  further  along  than  others.  This  motivates 
the  introduction  of  “local  clocks”  since  there  is  no  true  notion  of  global  time  in 
an  asynchronous  distributed  simulation.  When  a  process  sends  a  message,  it 
appends  to  the  message  a  time  stamp  which  bears  the  value  of  its  local  clock  at 
the  instant  when  the  message  was  generated.  The  contents  of  this  message  may 
alter  the  evolution  of  the  receiving  process;  however  the  time-stamp  will  gener¬ 
ally  not  match  the  receiving  process’s  local  clock.  This  means  that  the  message 
may  be  in  the  “past”  or  “future”.  Future  events  pose  little  problem,  however 
handling  past  events  requires  that  the  process  must  undo  or  “roll  back”  its  sim¬ 
ulation  to  a  point  where  its  local  time  clock  is  earlier  than  the  time  at  which  the 
message  was  generated.  This  necessarily  implies  that  the  process  must  store  a 
history  of  its  evolution  so  that  it  may  always  roll  the  simulation  back.  Roll  back 
becomes  expensive  when  many  processes  are  involved  because  the  rollback  of 
one  process  may  in  turn  trigger  the  rollback  of  another  and  so  on,  resulting  in 
what  is  called  a  cascading  rollback. 

An  active  area  of  research  in  DDES  concerns  finding  the  best  way  to  minimize 
rollback.  Approaches  fall  into  two  categories:  conservative  and  optimistic.  Con¬ 
servative  approaches  only  allow  each  process  to  proceed  up  until  a  certain  safe 
time  [Fujimoto  2000]  (a  time  before  which  it  can  be  guaranteed  no  rollback  will 
be  necessary).  Some  processes  must  wait  idle  upon  reaching  this  point  in  their 
simulation  until  other  processes  catch  up  resulting  in  underutilized  resources. 

In  Jefferson  [1985],  the  Time  Warp  Algorithm  was  introduced  which  was  the 
first  optimistic  simulation  methodology.  In  the  optimistic  approach  each  pro¬ 
cess  is  allowed  to  advance  at  its  own  rate  hence  no  resources  sit  idle;  however 
frequent  rollback  is  possible.  Jefferson’s  algorithm  included  an  elegant  mech¬ 
anism  called  antimessaging  to  reliably  deal  with  rollback.  Since  then  there 
has  been  an  explosion  of  research  in  the  field  exploring  the  tradeoff  between 
idle  resources  and  the  cost  of  rollback  (see,  e.g.,  Lubachevsky  [1989],  Steinman 
[1992],  or  Carothers  et  al.  [1999]).  Many  are  still  built  on  the  Timewarp  algo¬ 
rithm.  They  explore  such  issues  as  the  best  strategy  for  periodically  exchanging 
information  to  minimize  rollback.  Many  also  consider  topics  specific  to  dis¬ 
tributed  computation  such  as  network  delays  and  bandwidth.  It  is  obvious 
that  this  work  bears  relation  to  multi-agent  hybrid  systems;  however  since 
the  systems  are  purely  discrete  there  are  no  challenges  related  to  integrating 
the  differential  equations.  Also,  the  challenges  associated  with  implementing 
a  simulation  algorithm  on  multiple  processors  are  beyond  the  scope  of  this 
article. 
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Mirtich  [2000]  introduced  a  multibody  Newtonian  mechanics  simulator 
for  graphics  applications  which  is  built  upon  the  Time  Warp  algorithm 
[Jefferson  1985].  In  this  framework,  each  rigid  body  is  considered  a  process 
with  its  own  local  clock;  collisions  and  other  changes  in  contact  states  (rolling, 
sliding,  etc.)  play  the  role  of  messages  or  discrete  events.  This  framework  per¬ 
mits  the  use  of  an  asynchronous  integration  procedure.  Although  there  is  some 
additional  overhead  associated  with  the  method  as  the  number  of  bodies  in¬ 
creases,  the  performance  of  the  method  appears  to  be  superior  to  traditional 
synchronous  rigid  body  simulators  by  approximately  one  order  of  magnitude. 
It  was  shown  that  these  savings  increase  as  the  number  of  agents  grow.  The 
goal  of  this  work  is  identical  in  spirit  to  Mirtich’s — reduce  the  computational 
effort  required  to  simulate  large  multi-agent  systems  through  asynchronous 
simulation  and  careful  handling  of  discrete  events.  The  key  difference  is  that 
Mirtich’s  work  exploits  many  unique  features  of  Newtonian  Mechanics  and  uses 
collision  detection  algorithms,  specific  to  polyhedral  bodies,  to  detect  “events”. 
This  specificity  allows  the  simulator  to  achieve  only  slightly  subrealtime  perfor¬ 
mance;  however,  it  is  not  general  enough  to  apply  to  arbitrary  hybrid  systems. 
Similar  issues  are  addressed  in  Nicol  and  Perrone  [2000]  which  addresses  the 
special  structure  of  wireless  communication  systems. 

A  second  important  related  work  is  Hur  and  Lee  [2002],  which  introduced 
an  algorithm  for  distributed  (and  asynchronous)  multi-agent  hybrid  system 
simulation.  Several  heuristics  for  determining  the  frequency  with  which  global 
state  information  should  be  broadcasted.  The  most  obvious  difference  between 
this  work  and  ours  is  that  the  algorithms  and  implementation  are  distributed 
in  Hur  and  Lee  [2002],  and  therefore  different  issues  must  be  addressed.  A 
more  subtle,  yet  conceptually  important,  distinction  is  that  the  algorithms  pre¬ 
sented  in  Hur  and  Lee  [2002]  represent  optimistic  approaches  to  the  simulation 
problem — meaning  that  they  accept  possibly  frequent  rollback  in  the  local  time 
clocks.  Our  work  presented  here  can  be  viewed  as,  to  borrow  from  the  DDES  lit¬ 
erature,  a  conservative  approach  to  the  problem.  We  feel  that,  unlike  DDES,  the 
only  way  to  truly  guarantee  that  no  events  are  missed  in  hybrid  simulation  is  to 
only  permit  the  local  clocks  to  advance  independently  to  the  extent  that  it  can  be 
shown  that  no  events  will  be  missed.  In  many  ways  our  work  and  Hur  and  Lee 
[2002]  can  be  viewed  as  complementary  approaches  to  multi-agent  simulation. 

4.  APPROACH 

In  this  section,  we  describe  the  idea  behind  our  approach  to  selecting  step  sizes 
for  the  asynchronous  integration  in  such  a  way  as  to  reliably  detect  events, 
with  minimal,  if  any,  rollback.  We  defer  details  of  the  implementation,  as  well 
as  other  considerations  for  selecting  step  sizes  until  Section  5. 

Given  the  Multi-Agent  Hybrid  System 

f  =  fix1),  i  =  l,...,N  (6) 

gix1,  . . .  ,xN)  <  0,  (7) 

where  superscripts  index  the  N  agents.  Let  xl  e  R"'  be  the  state  of  Agent  i, 
fl  :  Rn‘  — »•  Rn‘  the  time  derivative  of  state  xl ,  and  g  :  Rni  x  •  •  •  x  R"N  -*  R  is 
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the  event  function.  Assume  that  f'  is  continuous  as  long  as  gix1, . . . ,  xN)  <  0 
and  that  g  itself  is  smooth.  Note  that  with  no  loss  of  generality,  one  may  assume 
that  g (x 1 , . . . ,  xN)  is  linear  in  the  states  since  systems  with  a  nonlinear  event 
function  can  be  converted  to  an  equivalent  system  with  a  linear  event  function 
by  appending  an  extra  state  variable  as  follows  [Gear  and  Osterby  1984]: 


xl 

—  fl{xl)  i  —  1, . 

..,N 

(8) 

xN+1 

■  fNf 

(9) 

g 

=  xN+1  <  0, 

GO) 

where  X  is  the  global  state  vector.  For  now,  we  assume  all  event  functions  are 
linear. 

Since  differential  equations  are  numerically  integrated  using  difference 
equations  and  the  states  are  approximated  at  a  discrete  set  of  points  in  time 
T  e  {t\,  t<i, . . . ,  tk,  ■  . .},  let  subscripts  denote  the  step,  or  iteration,  number.  Asyn¬ 
chronous  simulation  implies  there  is  no  single  global  value  of  time,  tk ,  instead 
the  simulation  of  each  agent  proceeds  independently,  whenever  possible.  As 
such,  local  times  and  time  steps,  must  be  defined  for  each  agent  and  they 
are  also  denoted  using  superscripts.  Therefore,  any  quantity  with  a  super¬ 
script  i  and  subscript  k,  refers  to  agent  i  and  has  been  evaluated  at  tlk.  For 
example,  xj,  =  xl(t'k)  and  f'k  —  f!(xl(t'k)).  Also  define  the  integration  time  step 
l,i  —  fi  _  /' 
ak  —  Lk  Lk-V 

As  in  Esposito  et  al.  [2001b],  we  choose  to  use  Predictor — Corrector  numerical 
methods  to  integrate  the  system  (see  any  standard  monograph  on  numerical  in¬ 
tegration  such  as  Ascher  and  Petzold  [1998]).  Applying  the  Predictor  difference 
equation  to  the  state 

m 

xl+i  =  4  +K+iJ2Pj(hk+i)fk-j+i>  <n> 

7  =  1 

where  the  /fs  are  weighting  functions  which  are  polynomial  functions  of  the 
step  size. 

When  Eq.  (11)  is  substituted  into  g(x )  the  change  in  the  value  of  the  event 
function  with  each  iteration  is  given  by, 

(m 

4 + hl+i  fr  (/i*+i)  fk-j+i>  ■  ■  ■ 

7  =  1 

m  \ 

xk  +  K+ 1  £  4  (hk+i)  fk-j+i )  (12) 

or,  by  linearity, 

dg  m 

Sk+i  gk  +  4+1  Pj  (hk+i)  fk-j+i 

7=1 

r,  m 

+  /l*+1^2  11,  Pj  ( hk+l )  fk-j+1  4 — 

7  =  1 
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rj  m 

+  hk+l  E  Pj  (hk+l)  fk-j+1  ■  (13) 

j= 1 

Using  a  control  system  analogy,  Eq.  (13)  is  analogous  to  a  system  with  a 
single  output,  g,  and  N  inputs,  h\  v  ...,h^v  Esposito  et  al.  [2001a].  Now 
our  goal  is  to  select  these  N  inputs  in  such  a  way  as  to  properly  “steer”  the 
system  to  any  zeros  of  g.  However,  it  is  important  to  point  out  that  simply 
driving  g(xA(tA), . . . ,  xN(tN))  -»  0  is  not  sufficient  to  detect  an  event.  In  fact 
g  —  0  does  not  correspond  to  a  physical  event  unless  th  1  =  t2+1  —  •  •  •  =  tj*v  If 
g  —  0,  but  tlk+1  /  tJk  v  this  implies  that  the  two  agents  have  passed  though  the 
same  point  on  the  event  surface  but  at  different  local  times.  Such  a  situation 
is  not  a  state  event,  it  is  simply  an  artifact  of  the  asynchronous  simulation. 
Returning  to  the  automated  vehicle  example,  if  the  guard  signified  collisions 
between  the  vehicles  (i.e.,  g(x)  is  separation  distance),  the  case  in  which  g  —  0 
but  tl+l  /  t'f1  would  imply  the  cars  passed  through  the  same  point  in  the 
configuration  space  at  different  times,  which  obviously  does  not  constitute  a 
collision. 

Fortunately,  the  fact  that  there  are  more  inputs  than  outputs  implies  that 
there  is  some  freedom  in  the  selection  of  the  step  sizes.  This  latitude  may  be 
used  to  accommodate  secondary  criteria.  Thus  N  —  1  additional  inequalities  are 
added  to  the  original  event  criteria,  which  we  call  synchronization  constraints. 
Since  the  labelling  of  the  agents  is  arbitrary,  assume  that  tji+1  >  t2+1  >  •  •  •  > 
t^+v  The  N  —  1  synchronization  functions  would  be  defined  as 

Tk+i  =  4+1  ~  4+1  <  °>  i  =  2,...,N,  (14) 

each  measuring  the  extent  to  which  an  individual  agent’s  time  clock  lags  the 
leading  agent’s  clock. 

Physically  meaningful  events  are  now  defined  as 

ig  =  0)  A  (r|+1  =  0)  A  ■  ■  ■  A  (t^+1  =  0).  (15) 

One  must  then  select  h\+1, . . . ,  h£+1  so  that  the  simulation  will  precisely  land 
on  the  point  in  time  when  all  of  these  conditions  are  true.  Collectively,  the 
dynamics  of  the  event  and  synchronization  functions  can  be  written  as  follows. 


"  g  " 

"  g  " 

9 

9 

r 

X 

_rN  _ 

k+i 

_rN  _ 

(16) 


^j:fN^hN)i 

-  h 1  - 

-1  1 

0.. 

0 

h2 

-1  0 

i.. 

0 

-1  o 

0.. 

_hN  _ 

1 

k 

-ik+i 
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Inspired  by  the  feedback  control  analogy  our  step  size  selection  scheme  is  es¬ 
sentially  a  discrete  time  version  of  the  feedback  linearization  technique  from 
nonlinear  control  theory,  which  uses  a  judicious  selection  of  the  input  to  cancel 
the  system’s  nonlinear  dynamics,  rendering  it  a  stable  linear  system. 

Theorem  4.1.  By  selecting  choosing  a  value  of  the  “gain”  0  <  y  <  1,  defining 

K+i  =  (y~  r>Tk  +  h\+ 1>  i  =  2,...,N  (17) 

and  selecting  h\+l  as  the  smallest  positive  real  root  of  the  following  polynomial 

o  m 

(1  ~Y)gk  +  hl+1^J2Pj(hl+i)fk-j+i 
l=i 

r)  a  171 

+  ((y  -  Dt*  +  h\+i)  ^  Z  Pj  (( Y  -  Dt*  +  h\+ 1)  A-j+1  (18) 

7=1 

+  •••  +  ((/-  l)rf  +  /4+1)  ^  Z  £/  ((y  -  l)rf  +  ^+1) 

1=1 

all  events  will  be  properly  detected  ( within  the  tolerance  level  of  the  integration) 
without  having  to  roll  the  simulation  back  ( i.e.,  decrement  k). 

Proof.  Substituting  Eq.  (17)  into  Eq.  (16)  produces  the  following  difference 
equations  for  the  synchronization  function 

gk+i  =  gk  +  h\+ 1-£I  J2  Pj  (hl+ 1)  fk-j+i  (19) 

7=1 

+  ((y  -  Dt*  +  ^*+i)  Z  Pj  (( y  -  1)tI  +  /l*+i)  A-i+i  +  •  ■  • 

i=i 

o  m 

+  ((y  -  Drf  +  fyt+i)  ^  Z  Pj  _  1)rf  +  /l<U)  /ft-i+i  (20) 

1=1 

T*+i  =  Y?lk,  i  —  2, . . . ,  N .  (21) 

Further  substituting  the  roots  of  Polynomial  (18)  for  h\+l  will  yield 

gu+ i  =  ygk  (22) 

4+i  =  Y  4-  (23) 

The  solution  to  the  above  difference  equations  is 

gk  =  (y)Vo  (24) 

T!  =  (Y)k  4  (25) 

Provided  0  <  y  <  1,  as  k  -»  oo,  gk  ->  0  and  r|  ->  0  for  all  i  =  2, ... ,  N, 
implying  that  the  simulation  will  terminate  at  the  event.  Furthermore,  if  ini¬ 
tially  g,  r1, . . . ,  rN  <  0,  this  equilibrium  point  is  approached  strictly  from  the 
right.  This  implies  that  there  is  no  danger  of  advancing  the  local  clocks  too  far, 
eliminating  the  need  for  rollback. 
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Several  observations  are  worth  pointing  out. 

(1)  The  step  size  selection  suggested  above  represents  the  ideal  step  size  based 
on  event  considerations.  In  practice,  the  step  size  will  occasionally  be  limited 
further  to  reduce  truncation  error  in  the  integration. 

(2)  Polynomial  (18)  is  essentially  the  most  accurate  extrapolation  polynomial 
that  can  be  constructed  which  is  consistent  with  the  underlying  integration 
accuracy.  Since  underlying  state  estimates  are  generated  by  a  forward  in¬ 
tegration  scheme  it  does  not  make  senses  to  use  the  negative  real  root(s) 
of  this  polynomial.  Obviously,  complex  roots  have  no  physical  meaning  as 
time  steps. 

(3)  If  there  are  no  positive  real  roots  of  the  polynomial,  then  from  the  point  of 
view  of  event  detection,  the  acceptable  step  size  is  unbounded. 

(4)  If  y  =  0,  the  system  will  converge  in  a  single  iteration. 

(5)  If  the  event  function  is  nonlinear  the  technique  is  still  applicable,  however 
the  statement  regarding  the  need  for  rollback  must  be  relaxed.  The  above 
choice  of  step  size  will  still  cancel  the  first  and  most  dominant  terms  of 
the  Taylor  Series  expansion  of  Eq.  (12),  however  the  uncancelled  higher- 
order  terms  may  act  as  forcing  functions  and  introduce  the  possibility  of 
overshooting  the  exact  event  time  by  at  most  one  step.  This  is  discussed 
further  in  Esposito  [2002]. 


5.  SIMULATION  ALGORITHM 

The  multi-agent  simulation  algorithm,  shown  below,  is  referred  to  as  MAsim. 
Figure  1  illustrates  the  simulation  algorithm  graphically;  while  the  pseudo-code 
in  Algorithms  1  and  2  provide  more  detail.  The  process  illustrated  in  Figure  1 
resembles  a  traditional  hybrid  system  simulator,  using  a  Predictor-Corrector 
integration  algorithm,  in  many  ways.  First  an  initialization  phase  is  required 
because  Predictor-Corrector  integrators  require  a  history  of  derivative  values. 
Different  implementations  achieve  this  start  up  process  in  different  ways.  Note 
that  the  user  specifies  e,  the  integration  error  tolerance  which  should  not  be 
exceeded.  The  step  size  computation,  as  discussed  in  Section  4  is  done  by  the 
routine  Algorithm  2,  called  MAstepSelect. 

An  important  issue  in  asynchronous  simulation  is  deciding  the  order  in  which 
the  individual  agents  should  be  advanced.  In  this  case,  at  each  iteration  in  the 
simulation  the  agent  whose  local  time  clock  lags  the  most  is  always  selected  to 
be  integrated.  In  the  case  of  a  tie,  an  agent  may  be  chosen  at  random  among  the 
agents  with  the  smallest  local  times.  At  the  end  of  each  iteration,  the  agents 
are  relabelled  so  that  $+1  ^  *£+1  >■■■>$+!  to  ensure  that  r^+1  <  0.  Therefore, 
at  each  iteration,  agent-N  is  always  selected  for  advancement.  This  process 
of  consistently  selecting  the  agent  with  the  lagging  time  clock  has  the  effect 
of  bounding  the  maximum  amount  of  time  by  which  the  local  time  clocks  may 
disagree.  In  effect, 

max  [rl]  <  hmax,  (26) 

i.=2,...,N 
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Fig.  1.  The  overall  simulation  algorithm  MAsim  presented  in  Section  5. 


where  hmax  is  the  prespecified  maximum  allowable  step  size.  Since  the  event 
detection  scheme  is  essentially  based  on  solving  for  the  roots  of  an  extrapolation 
polynomial  at  each  step,  bounding  the  extrapolation  interval  helps  to  ensure 
meaningful  estimates  of  the  event  time. 

In  Figure  1,  everything  inside  of  the  dashed  box  is  accomplished  using  the 
standard  integration  technique  along  with  the  Predictor-Corrector  method. 
After  each  new  state  is  constructed,  the  approximation  error  in  the  result  is 
estimated  and  compared  to  a  threshold.  This  information  is  used  to  decide  if 
the  result  should  be  accepted  or  rejected  as  well  as  to  determine  what  herr 
should  be  to  keep  the  integration  error  within  the  acceptable  range  for  the  next 
iteration. 

Finally,  it  must  be  determined  if  an  event  has  occurred.  This  determination 
must  be  made  using  our  new  requirement  for  event  occurrence,  namely  — 
0  /\  r|  =  0  A  •  •  •  A  T}?  —  0.  If  indeed  one  has  occurred,  the  simulation  is  stopped, 
all  mode  switches  are  enabled  and  the  simulation  restarts  in  the  new  mode. 
Note  that  the  initialization  process  must  be  repeated  since  the  definitions  of 
the  derivatives  have  changed.  If  instead,  it  turns  out  that  no  event  occurred, 
the  agents  are  relabelled  so  that  tk  >  t2  >  •  •  •  >  t?  to  ensure  that  4  <  as 
discussed  above. 
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Algorithm  1  MAsim :  Multi-agent  simulation  algorithm 


0.  Given  initial  conditions  a:1(0), . . .  x'v(0);  user  defined  values  for  A,,,,,,,  hmax,  and  e;  and 
initial  estimates  for  h\rr  . . .  h^rr. 

1.  Initialize  the  integration  algorithm;  let  k  =  1. 

2.  Compute  h %  =  MAstepSelectC//, . . . ,  f^_m,  ■  ■  ■ 


fN  h  h  -  hN  —) 

/  k-mi  ^maxj  ,Lerr ’  a* 


A" 


3.  Integrate  Agent  Ar  using  /i^  via  the  Predictor  Corrector  Algorithm.  Estimate  error  e‘ 

if  eAf  <  e  then 

Store  result. 

Compute  hgrr  for  next  iteration. 


k  <-  k  +  1. 

else 

Reduce  h £ . 
Goto  3. 

end  if 


if  g  =  0  a  t2  A  ■  ■  ■  A  zN  then 
Switch  modes. 

Goto  1. 
else 

Relabel  agents  1, . . . ,  N  in  order  of  descending  values  of  t'  . 
Goto  2. 

end  if 


Algorithm  2  MAstepSelect:  The  step  size  selection  subroutine  for  MAsim. 


hN  =  MAstepSelect!  , . . . ,  fl_m, hmax,  hmin,  h*r,  Vg )  { 

for  i=  1:N  do 
T 1  =tl  -  t1. 

end  for 

Z  =  Roots  of  Polynomial  (18). 

R  =  {r|r  e  Z ,  Im(r )  =  0,  Re{r)  >  0). 

if  R  —  JD  then 

hg  =  oo. 

else 

hi  =  min(R). 
hg  =(y  -  1  )zjy  +  h1g. 

end  if 

hN  =  max[/imin,  minftmaY,  min[/ij,  h^rr\)].  ) 


Turning  our  attention  to  the  subroutine  that  selects  the  integration  step 
sizes,  there  are  several  issues  warranting  explanation.  Note  that  even  though 
routine  seeks  to  compute  hN ,  h1  must  be  computed  first  because  hN  is  defined 
in  terms  of  hN  as  in  Eq.  (17).  Step  sizes  with  a  subscript  err  denote  candidate 
step  sizes  computed  on  the  basis  of  controlling  the  truncation  error;  while  the 
subscript  g  denotes  a  step  size  computed  based  on  event  detection  criteria.  As 
discussed  in  Section  4,  the  ideal  step  size  from  the  point  of  view  of  event  detec¬ 
tion  is  computed  according  the  smallest  real  positive  roots  of  the  polynomial.  If 
there  are  no  such  root,  then  hg  is  set  to  infinity.  This  implies  that,  from  an  event 
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Fig.  2.  Two  automated  vehicles  driving  in  spiral  trajectories. 


detection  perspective,  there  is  no  constraint  on  the  step  size.  Finally,  note  that, 
as  with  any  simulation,  event  proximity  is  not  the  only  consideration  when 
choosing  a  step  size.  Typically,  there  are  user  defined  minimum  and  maximum 
step  sizes,  hm in  and  /imax,  as  well  as  an  estimate  of  the  ideal  step  size  based  on 
the  truncation  error  herr.  The  correct  choice  is 

hk+i  =  max |  /?.  m j n ,  min(/imax,  mini  h^rr ,  hg  I )  I .  (27) 


6.  EXAMPLES 

In  this  section,  the  effectiveness  of  the  algorithm  is  demonstrated  on  the  two 
car  example  used  to  motivate  the  problem  in  Section  1,  also  illustrated  in 
Figure  2.  Events  occur  when  the  two  autonomous  vehicles  come  within  R 
feet  of  each  other,  in  practice  this  condition  would  signal  an  imminent  col¬ 
lision  and  a  mode  switch  would  activate  an  emergency  collision  avoidance 
controller. 

The  dynamics  for  agent  i  =  1,2  are  given  in  Eq.  (2)  and  Eq.  (3);  while  the 
guard  is  given  in  Eq.  (5).  The  vl  and  o/  are  state  feedback  laws  (i.e.,  func¬ 
tions  of  xl,  yl  and  0l)  and  are  selected  in  such  a  way  as  to  steer  the  vehi¬ 
cles  in  spiral  trajectories.  However,  Car  1  travels  much  faster  than  Car  2. 
This  example  was  simulated  using  the  algorithm  presented  earlier  for  two 
scenarios-each  with  a  different  value  for  the  of  initial  states  and  separation 
distance. 

In  the  first  example  (see  Figure  3  and  4)  the  value  of  D  is  small  enough  that 
the  robots  eventually  collide.  The  path  of  the  robots  in  the  plane  is  shown 
in  Figure  3  (left).  Figure  3  (right)  is  a  plot  of  the  value  of  the  event  func¬ 
tion  versus  the  step  number  k.  The  value  peaks  first  around  step  70,  when 
Robot  1  completes  a  half  of  a  rotation;  finally  at  step  number  107  the  event 
detection  criterion  becomes  active  and  the  step  sizes  are  selected  in  such 
a  way  that  g  ->  0  exponentially.  Further  insight  is  gained  by  examining 
Figure  4  (left),  which  shows  the  step  sizes  used  for  each  agent  as  a  func¬ 
tion  of  time.  Note  how  h 1  and  h2  vary  independently  until  t  ~  9.8  when  the 
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Fig.  3.  The  trajectory  of  the  two  robot  vehicles  in  the  plane  is  plotted  (left);  it  is  of  interest  to 
detect  when  they  come  within  R  feet  of  each  other.  The  value  of  the  event  function  (separation 
distance)  as  a  function  of  step  number  is  shown  (right).  The  event  function  value  converges  to  zero 
exponentially. 


Fig.  4.  Step  sizes  used  in  the  first  example  are  shown  (left),  h 1  and  h 2  are  selected  independently 
away  from  the  constraint  but  are  brought  into  synchronization  when  an  event  is  impending.  The 
value  of  r  (right),  which  is  a  measure  of  the  discrepancy  between  the  two  local  time  clocks  (i.e., 
asynchrony),  rapidly  decreases  in  magnitude  as  the  event  is  approached. 


event  detection  criteria  begins  constraining  the  step  sizes;  whereafter  they 
both  begin  adjusting  the  step  sizes  so  as  to  synchronize  the  two  local  clocks. 
This  is  further  illustrated  in  Figure  4  (right),  which  plots  the  history  of  the 
synchronization  function  r.  Its  value  rapidly  approaches  zero,  at  which  point 
t 1  =  t2 . 

In  the  second  case  (see  Figures  5  and  6),  the  value  of  the  initial  separation 
distance  is  large  enough  that  the  vehicles  do  not  collide,  but  they  do  come  very 
close  to  one  another.  The  trajectories  of  the  two  vehicles  are  shown  in  Figure  5 
(left).  There  is  a  near  miss  halfway  through  Agent-l’s  second  rotation.  The 
history  of  the  constraint  in  Figure  5  (right)  shows  that  the  value  of  g(x)  comes 
very  close  to  zero  around  step  number  150.  During  this  period,  the  value  of  r, 
Figure  6  (left),  decreases  in  preparation  for  a  possible  event.  Figure  6  (right) 
shows  how  the  step  sizes  were  decreased  to  slow  down  the  simulation.  However 
once  the  two  vehicles  pass  each  other  and  it  becomes  apparent  that  no  collision 
will  occur,  the  step  sizes  quickly  increase  and  the  two  local  clocks  are  allowed 
to  fall  out  of  synchronization  again.  Note  that  in  these  simulations,  y  —  0.5  so 
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Fig.  5.  The  trajectory  of  the  two  robot  vehicles  in  the  plane  for  the  second  example  (left).  The 
initial  condition  is  selected  such  that  the  two  vehicles  do  not  come  within  R  feet  of  one  another. 
However  they  do  approach  quite  closely.  The  value  of  the  event  function  as  function  of  the  step 
number  is  shown  (right). 


Fig.  6.  The  step  sizes  used  in  the  second  example  are  shown  (left).  Again,  h 1  and  h 2  are  selected 
to  synchronize  the  simulation  when  it  seems  an  event  may  occur.  Once  it  is  apparent  that  this  is 
not  the  case,  the  simulation  returns  to  asynchrony.  The  magnitude  of  r  (right),  therefore  decreases 
during  this  period  but  soon  increases  again. 


that  the  “slow  down”  and  synchronization,  as  the  state  approached  the  guard 
surface,  could  be  more  easily  illustrated.  In  practice  a  value  closer  to  zero  would 
produce  rapid  convergence. 

7.  EFFICIENCY  ANALYSIS 

Recall  that  the  primary  motivation  for  using  asynchronous  simulation  is  to  im¬ 
prove  the  efficiency  of  the  simulator.  In  this  section,  the  predicted  performance 
of  the  asynchronous  algorithm  is  modelled  statistically  under  some  assump¬ 
tions  and  compared  to  that  of  the  more  traditional  asynchronous  algorithm. 

7.1  The  Traditional  Algorithm 

For  the  sake  of  comparison,  the  traditional  synchronized  algorithm  is  reviewed 
here.  The  traditional  algorithm  makes  no  distinction  between  single  agent 
systems  and  multi-agent  system.  In  both  cases  it  treats  the  entire  state  vector 
as  if  it  belongs  to  a  single  system  by  using  a  single  global  step  size;  we  refer 
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Algorithm  3  SAsim:  The  traditional  (i.e.,  single  agent)  hybrid  system  simulation 
algorithm. 


Given  initial  conditions  x'(0), . . . ,  xN(0);  user  defined  values  for  hrmn,  hmax,  and  e;  and 
initial  estimates  for  herr. 

1.  Initialize  the  integration  algorithm.  Let  k  =  1. 

2.  Compute  hk  =  SAstepSelect]  f£, fl_m, . . . ,  , . . . ,  f£Lm,  /imax,  n,  h^r,  Vg). 

3.  fori  =  1:N  do 
Integrate  xl  using  hk . 

Estimate  error  e1. 

end  for 

e  =  maxte1, . . . ,  eN). 

if  e  <  e  then 

Store  result. 

Compute  herr  for  next  iteration. 

k  •«—  k  +  1. 

else 

Reduce  hk  . 

Goto  3. 

end  if 

if  g  =  0 

Switch  modes. 

Goto  1. 
else 
Goto  2. 
end  if 


Algorithm  4  SAstepSelect:  The  step  size  selection  algorithm  used  by  SAsim,  the 
traditional  simulation  algorithm 

h  =  SAstepSelect!  ,  ■■■,  fk-m,  ■  ■  ■ ,  fk,  ■  ■  ■ ,  fk-m,  hm*x,  h^,  herr,  Vg)  { 

Z  =  Roots  ( P(h )). 

R  =  {r\r  e  Z ,  Imir)  =  0,  Re{r )  >  0). 

if  R  =  JO  then 

hg  =  oo. 

else 

hg  =  min(r). 

end  if 

Select  final  step  size  h  =  maxl/imin,  min(/ig,  herr,  /imax))- 


to  this  as  the  Single  Agent  Algorithm  or  SAsim.  The  notation  follows  the  same 
convention  as  in  the  Multi-agent  Algorithm  MAsim. 

While  any  suitable  integration  method  or  extrapolation  polynomial  can  be 
used,  since  the  goal  of  this  discussion  is  to  illuminate  the  pros  and  cons  of 
asynchronous  simulation  with  respect  to  the  approach  of  using  a  single  global 
step  size  we  assume  the  numerical  integration  in  both  cases  is  performed  using 
the  same  algorithm  (Predictor-Corrector)  and  a  very  similar  approach  to  event 
detections  such  as  that  of  Esposito  et  al.  [2001b].  The  polynomial  P(h)  which 
appears  in  the  step-select  routine  within  the  SAsim  algorithm  is  essentially 
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a  “single  agent”  version  of  Eq.  (18) 


(1  -  y)gk  +  hk+ i 


dx1 


j+ 1 


j= i 


9g 


/*  ,  P.i (hk+i)  fk- 7+i 


f=i 


(28) 


Note  that  the  superscript  on  li  has  been  dropped  because  there  is  a  single  global 
step  size;  however  the  polynomial  has  the  same  effect  of  feedback  linearizing 
the  guard  dynamics. 

There  are  two  primary  differences  between  the  SAsim  and  MAsim  algo¬ 
rithms.  In  SAsim,  a  single  ideal  step  size  from  the  point  of  view  of  event  detec¬ 
tion  hg  is  computed  using,  while  in  MAsim  the  event  dynamics  in  Eqs.  (18)  and 
(17)  are  used  to  compute  N  distinct  step  sizes,  one  for  each  agent.  Obviously, 
doing  this  requires  significantly  more  overhead  in  computing  the  various  r’s, 
h’s  and  relabelling  the  agents.  The  underlying  polynomial  used  in  MAsim  is 
also  slightly  more  complex. 

The  second  difference  is  that  in  SAsim  a  single  step  size  from  the  point 
of  view  of  the  truncation  error  is  computed.  This  step  size  is  based  on  the 
maximum  of  all  the  agents  truncation  error.  This  results  in  using  the  minimum 
step  size  across  all  agents.  Obviously,  this  implies  the  SAsim  algorithm  will,  in 
general,  use  smaller  time  steps  and  therefore  requires  more  steps  to  complete 
the  simulation. 


7.2  Average  Complexity  Analysis 


Overall  MAsim  requires  more  computation  per  iteration;  while  SAsim  requires 
more  iterations  to  complete  a  simulation.  In  this  section,  we  attempt  to  quantify 
this  tradeoff  by  studying  their  affect  on  the  quantity,  which  we  call  the  Speed-up 
Factor,  defined  as 


S  = 


Cma 
CsA  ’ 


(29) 


where  Cma  is  the  computational  cost  (floating  point  operations  or  CPU  time) 
required  to  simulate  a  specific  system  using  MAsim  and  C$a  is  the  cost  of 
simulating  the  same  scenario  using  SAsi?n,  assuming  all  parameters  are  equal, 
such  as  the  integration  tolerance.  Values  of  S  significantly  smaller  than  one 
imply  that  a  great  speedup  is  obtained  using  MAsim. 

Obviously,  S  is  entirely  problem  (and  implementation)  dependent  therefore 
we  choose  to  instead  study  the  average  Speed-up  Factor.  The  assumption  is  that 
since  many  complex  factors,  which  cannot  all  be  modelled,  influence  the  step 
size  selection  criteria  we  choose  to  make  the  reasonable  modelling  assumption 
that  the  step  size  that  the  algorithm  selects  is  a  uniformly  distributed  random 
variable  h’  e  |/irrun,  /imax  I-  This  assumption  enables  us  to  focus  on  the  general 
effects  of  the  tradeoff  between  having  a  higher  overhead  in  the  case  of  MAsim 
versus  being  forced  to  use  smaller  step  sizes  in  the  case  of  SAsim.  Note  that  the 
analysis  can  be  repeated  using  any  distribution  on  the  step  sizes;  however,  the 
uniform  distribution  is  chosen  because  it  does  not  clearly  favor  either  algorithm. 
In  that  sense,  it  is  arguably  the  most  “neutral”  distribution.  Further  discussion 
is  postponed  to  Section  8. 
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Since  Csa  =  E[hsA\CI  and  Cma  —  EUimaKCI  +  OH)  where  E[Jisa\  and 
E[Iima\  are  the  expected  values  of  the  step  sizes  chosen  by  the  two  algorithms; 
OH  is  the  additional  overhead  associated  with  the  MAsim  algorithm;  and  C/ 
is  the  cost  of  integrating  a  single  agent  through  a  single  time  step.  Note  that 
Cj  is  identical  in  both  algorithms,  provided  the  same  integration  technique  is 
used.  Therefore, 


_E[hSA]  f  Offl 

EUima]  1  Cj  ) 


(30) 


Savg  then  depends  on  two  factors:  the  ratio  of  the  expected  values  of  the  step 
size  selected  by  each  algorithm  E[hsA\/ E[Iima%  and  the  size  of  MA’s  overhead 
in  comparison  with  the  cost  of  integrating  though  a  single  step,  OH/Ci,  (re¬ 
ferred  to  as  the  relative  overhead ).  This  quantifies  the  discussion  in  Section  7.1. 


7.2.1  Expected  Values  of  Step  Sizes.  The  assumption  of  a  uniform  distri¬ 
bution  means  that  the  expected  value  of  the  step  size  for  the  ith  agent  is  simply 


Em 


dim  ax  h  min) 

2  +  Zlm  in 


(31) 


In  the  case  of  the  Multi-agent  Algorithm,  at  each  iteration  step  sizes  h1, ... ,  hN 
are  computed  and  each  agent  is  integrated  with  its  respective  step  size.  So 
EUima]  =  E[hl\  since  the  expected  value  does  not  depend  on  i.  The  distribution 
for  Iima  remains  the  same  regardless  of  the  number  of  agents. 

In  the  case  of  the  Single-agent  Algorithm  h1, . . . ,  hN  are  computed  in  the 
same  way;  however,  at  each  iteration,  a  single  global  time  step  must  be  selected 
for  all  agents  due  to  synchronization  requirements,  so 

Iisa  —  minf/i1, . . . ,  hN].  (32) 


Therefore,  one  would  expect  the  distribution  of  Iisa  to  be  biased  toward  small 
step  sizes  (see  Figure  7).  It  is  difficult  to  write  a  closed  for  expression  for  E[hsAi, 
however,  a  normalized  histogram  generated  numerically  appears  in  Figure  8. 
As  N  increases,  one  sees  that  the  frequency  with  which  small  step  sizes  are 
selected  increases  dramatically  in  the  Single-agent  Algorithm.  It  can  be  shown 
that  the  expected  value  of  this  distribution  takes  the  form 


E[IiSa\  »  Zimin  +  Cl  exp(— C21V), 

(33) 

where  Ci  and  C2  are  positive  constants.  As  a  result 

EUisa\  ~  Zimin  +  Cl  exp(— C2A0 

(34) 

E[hMA]  (^max  T-  ^min)/2 

Therefore, 

E[_hsA~]  2/imjn 

N  ^oo  E  [hMA]  ^max  “1“  ^-min 

(35) 

7.2.2  Computing  the  Relative  Overhead.  The  second  factor  influencing 
Savg  is  the  relative  overhead,  which  is  essentially  the  cost  of  completing  any 
extra  steps  in  MAsim  not  present  in  SAsim  such  as  computing  the  synchro¬ 
nization  events,  relabeling  the  agents,  and  calculating  extra  step  sizes.  The 
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Fig.  7.  The  expected  value  of  h  (normalized),  for  the  Single  Agent  Algorithm  as  a  function  of  the 
number  of  agents. 


Fig.  8.  A  histogram  depicting  the  frequency  of  selected  step  sizes  (normalized)  for  the  Single-agent 
algorithm,  for  several  values  of  N  (the  number  of  agents),  as  the  number  of  agents  becomes  large. 


complexity  of  each  these  components  of  the  overhead  is  O(N)  in  the  number  of 
agents, 

OH 

—  «  c3  +  C4N,  (36) 

where  C3  and  C4  are  positive  constants.  This  was  confirmed  experimentally. 
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Fig.  9.  The  iV-agent  “bumper  car”  example.  Events  occur  when  an  agent  collides  with  one  of  the 
walls  or  another  agent. 


Recall  that  the  relative  overhead  is  the  overhead  as  a  percentage  of  the  cost  of 
integrating  through  a  single  step  C/.  Naturally  C/  is  highly  problem  dependent 
but  it  is  at  least  10  floating  point  operations.  Problems  with  more  complicated 
expressions  for  f(x)  typically  have  higher  values  of  C/.  Such  problems  favor 
MAsim  even  further. 


7.2.3  Prediction  and  Discussion.  Using  (30),  (34),  and  (36)  on  can  then 
predict  the  behavior  of  Savg  as  a  function  of  the  number  of  agents, 


S 


avg 


(Amin  +  Cl  exp-c^XC3  +  C4N) 

(A  max  +  Amin)/2 


(37) 


where  the  particular  values  of  the  constants  are  problem  dependent. 


8.  NUMERICAL  EXPERIMENTS 

The  benchmark  problem  chosen  to  test  these  ideas  is  a  larger  version  of  the 
autonomous  vehicle  test  problem  first  discussed  in  Section  1  and  presented  in 
Section  6.  N  kinematic  car-like  agents  are  confined  to  a  rectangular  “arena”. 
Anytime  an  agent  comes  within  a  certain  distance  from  a  “wall”  or  another 
vehicle,  an  event  occurs  and  an  emergency  control  law  is  activated.  See  Figure  9. 
The  dynamics  are  given  by  Eq.  (2)  and  the  functions  v  and  a>  were  arbitrarily 

chosen.  When  a  collision  is  about  to  occur  with  either  another  agent  or  with 

the  boundary  of  the  arena,  the  control  law  switches  to  avoid  a  collision.  The 
“agent-wall”  collision  events  are 

g\  =  chv  -  xl  <  0  (38) 

g\  =  xl  -  crw  <  0  (39) 
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(40) 

(41) 


g\  =  yl  -ctw  <  o 

glb  =  cbw-  yl  <  0 

where  ciw  and  crw  and  are  the  x  coordinates  of  the  left  and  right  walls,  and  ctw 
and  cbw  and  are  the  y  coordinates  of  the  top  and  bottom  walls.  The  “agent- 
agent”  collision  event  is 

gij  =(xi-xj)2  +  (yi-yj)2-R2<  0,  Vj  /b  (42) 

where  R  is  the  range  at  which  the  emergency  control  law  kicks  in.  The  number 
of  differential  equations  is  3N  and  the  number  of  guards  is  AN  +  N  ■  ( N  —  l)/2. 

The  algorithms  SAsim  and  MAsim  were  coded  in  Matlab  so  as  to  be  as  similar 
as  possible  to  allow  a  fair  comparison.  Parameters  such  as  the  minimum  and 
maximum  step  size,  and  integration  tolerance  were  the  same  in  both  cases.  The 
actual  numerical  integration  algorithms  were  also  identical. 

Experiments  were  conducted  to  determine  how  the  speedup  factor  Savg 
varies  with  the  number  agents.  Recall  that  it  was  predicted  that  as  the  number 
of  agents  increases,  the  Single  Agent  algorithm  becomes  increasingly  likely  to 
select  smaller  step  sizes,  while  the  Multi-agent  Algorithm  is  not  biased  toward 
smaller  step  sizes.  Data  was  collected  for  systems  with  2  to  20  agents.  For  each 
number  of  agents,  the  system  was  simulated  for  10  randomly  generated  initial 
conditions,  using  a  relative  error  tolerence  of  10”3,  for  a  fixed  duration  using 
both  SAsi?n  and  MAsim.  The  number  of  floating  point  operations  used  by  each 
algorithm  for  a  given  scenario  was  measured  using  the  flops  command  in  Mat- 
lab.  The  raw  data  is  shown  in  Figure  10.  Note  for  most  of  the  190  test  scenarios, 
with  the  exception  of  2  of  the  two  agent  scenarios,  MA  outperformed  SAsim. 
In  the  case  when  N  —  20,  SAsim  required,  on  average,  double  the  amount  of 
computations  required  by  MA. 

A  second  factor  whose  impact  on  performance  was  studied  is  the  user- 
specified  integration  error  tolerance.  The  performance  improvement  associated 
with  MA  would  be  expected  to  decline  at  tighter  error  tolerances.  This  is  an¬ 
ticipated  because  tighter  error  tolerances  have  the  effect  of  biasing  the  step 
size  selection  mechanism  toward  smaller  step  sizes.  Each  of  the  190  scenarios 
mentioned  above  was  simulated  with  three  different  integration  errors  toler¬ 
ances,  e  =  10“3,  10"4  and  10-5,  for  a  total  of  570  experiments.  The  number 
of  floating  point  operations  used  by  each  algorithm  for  a  given  scenario  was 
measured. 

Based  on  this  data  the  speedup  ratio  for  each  scenario  was  calculated.  Savg 
for  a  given  number  of  agents  and  a  given  tolerance  was  computed  as  the  mean 
of  the  speedup  ratios  of  the  10  random  scenarios.  The  Savg  as  a  function  of  the 
number  of  agents  is  plotted  in  Figure  11.  Each  curve  displays  the  results  for 
a  different  value  of  the  integration  tolerance.  Tighter  integration  tolerances 
have  the  effect  of  biasing  the  algorithms  toward  selecting  smaller  step  sizes. 
The  dashed  lines  in  the  figure  shows  the  overall  trend,  qualitatively  predicted 
by  Eq.  (37). 

Clearly,  the  assumption  that  the  selected  step  size  is  a  uniformly  distributed 
random  variable  is  an  idealization  and  cannot  be  used  to  make  quantitative 
predictions  on  the  performance.  However,  the  data  in  Figure  11  and  the  shape 
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Fig.  10.  Sample  experimental  results  on  the  complexity  of  each  algorithm  as  a  function  of  the 
number  of  agents.  Squares  represent  SAsim  data  and  crosses  are  MAsim  data. 


Fig.  11.  Experimental  results:  Computational  cost  as  a  function  of  the  number  of  agents,  for 
various  integration  tolerances.  Dash  lines  represent  analytical  predictions;  while  the  solid  line  is 
the  geometric  mean  of  the  experimental  data. 


of  the  dashed  lines  obtained  by  fitting  Eq.  (37)  display  markedly  similar  quan¬ 
titative  trends.  The  stochastic  modelling  exercise  clearly  captures  the  essence 
of  the  tradeoff  at  work  between  step  size  selection  and  overhead.  It  also  indi¬ 
cates  that  a  significant  efficiency  gain  can  be  realized  in  practice  (in  the  case  of 
20  agents  a  factor  of  5  or  more). 
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Of  course,  the  actual  step  size  “distribution”  is  nonuniform,  which  warrants 
a  discussion  how  this  assumption  affects  the  analysis.  Distributions  which  are 
non-uniform  and  heavily  weighted  around  a  mean  will  certainly  bias  the  pre¬ 
diction  depending  on  how  the  mean  varies  across  all  agents.  If  all  the  agents 
have  the  same  distribution,  with  the  same  mean,  this  implies  that  the  agents 
are  all  varying  at  nearly  the  same  time  scale,  in  which  case  little  advantage  is  to 
be  gained  from  an  asynchronous  simulation.  On  the  other  hand,  if  all  agents 
have  the  same  general  distribution  but  their  means  do  not  coincide,  the  asyn¬ 
chronous  simulation  will  be  highly  favored  since  this  implies  that  some  agents 
are  consistently  “slow”  and  others  consistently  “fast”.  In  light  of  this,  a  uniform 
distribution  was  selected  since  it  is  does  not  necessarily  favor  either  of  these 
two  extremes,  representing  a  somewhat  neutral  assumption. 

9.  CONCLUSION 

A  major  cause  of  performance  loss  when  simulating  large  scale  systems  of  ODEs 
can  be  traced  back  to  the  requirement  that  a  single  global  step  size  is  tradi¬ 
tionally  used  despite  the  fact  that  the  various  individual  differential  equations 
in  the  system  may  each  warrant  a  very  different  choice  of  step  size.  The  only 
proper  choice  of  a  global  time  step  among  each  individual  equation’s  candidate 
time  step  is,  unfortunately,  the  minimum  of  these  candidates.  As  a  result,  the 
simulation  often  must  proceed  unnecessarily  slowly  using  the  smallest  possible 
step  size. 

In  situations  in  which  the  coupling  between  the  right-hand  sides  of  the  ODEs 
is  strong  this  cannot  be  avoided.  However,  we  consider  the  special  but  important 
case  of  decoupled  ODE’s  where  the  definition  of  the  derivatives  can  change,  or 
switch,  when  some  function  of  the  state  variables  changes  sign.  The  switching 
conditions  may  introduce  coupling  between  the  subsystems.  We  refer  to  such 
systems  as  multi-agent  hybrid  systems.  Examples  include  groups  of  automated 
highway  vehicles,  unmanned  aerial  vehicles,  mobile  robots,  etc. — systems  that 
are  nominally  physically  decoupled,  however,  through  sensing,  software  or  ac¬ 
tual  physical  collisions,  interactions  may  arise  which  precipitate  mode  changes 
within  the  system. 

In  this  article,  a  technique  is  introduced  for  simulating  multi-agent  hybrid 
systems  in  an  asynchronous  fashion,  that  is  without  resorting  to  a  single  global 
steps  size.  In  this  asynchronous  simulation  algorithm,  each  agent  has  its  own 
local  time  clock.  During  parts  of  the  execution  in  which  no  switches  occur,  the 
individual  time  clocks  are  allowed  to  advance  at  different  rates  as  dictated  by 
the  agent’s  dynamics,  allowing  the  maximum  utilization  of  the  available  com¬ 
putational  resources.  This  complicates  the  event  detection  problem  because  it 
is  no  longer  sufficient  to  simply  check  the  sign  of  the  guard  function  since  the 
individual  states  are  being  reported  at  different  points  in  time.  We  introduced 
and  algorithm  that  selects  N  integration  step  sizes  for  the  N  agents  in  such  a 
way  that  the  individual  time  clocks  will  synchronize  on-the-fly  when  some  rel¬ 
evant  subset  of  the  states  approach  a  guard  surface.  The  problem  is  analogous 
to  a  multi-input  multi-output  control  system,  for  which  the  inputs  are  the  N 
step  sizes  and  the  output  functions  are  the  guard  function,  along  with  N  —  1 
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synchronization  functions,  which  measure  discrepancies  between  the  local  time 
clocks. 

We  show  that  while  a  single  iteration  of  the  new  algorithm  is  somewhat  more 
expensive  due  to  the  overhead  associated  with  adaptive  synchronization,  not  us¬ 
ing  a  global  time  step  requires  fewer  total  iterations.  Because  the  performance 
is  somewhat  problem  dependent,  we  construct  a  probabilistic  model  of  the  step 
size  selection  scheme,  assuming  the  step  size  for  each  agent  is  a  random  vari¬ 
able.  Under  the  assumption  that  the  distribution  is  uniform,  it  is  shown  that 
average  complexity  of  the  multi-agent  algorithm  is  approximately  quadratic  in 
the  number  of  agents  being  simulated.  In  contrast,  it  is  demonstrated  through 
experiments  that,  as  the  number  of  agents  increases,  the  average  computation 
time  for  the  traditional  algorithm  seems  to  exponentially  approach  an  upper 
bound. 

Experiments  were  conducted  in  which  both  the  traditional  algorithm  and  the 
new  multi-agent  algorithm  were  used  to  simulate  a  problem  many  times  with  a 
large  number  of  randomly  generated  sets  of  initial  conditions.  Two  parameters 
were  varied:  the  total  number  of  agents  and  the  integration  error  tolerance. 
It  was  shown  that  the  predicted  relationship  between  the  number  of  agents 
and  the  computational  cost  savings  was  valid  and  that  requiring  finer  simula¬ 
tion  accuracy  serves  to  amplify  the  potential  savings  associated  with  the  new 
algorithm. 
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